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ABSTRACT 

In this paper we show how the Smoothed Particle Hydrodynamics (SPH) equations for ideal 
magnetohydrodynamics (MHD) can be written in conservation form with the positivity of the 
dissipation guaranteed. We call the resulting algorithm Smoothed Particle Magnetohydrody- 
namics (SPMHD). The equations appear to be accurate, robust and easy to apply and do not 
suffer from the instabilities known to exist previously in formulations of the SPMHD equa- 
tions. In addition we formulate our MHD equations such that errors associated with non-zero 
divergence of the magnetic field are naturally propagated by the flow and should therefore 
remain small. 

In this and a companion paper (Price & Monaghan 2003) we present a wide range of 
numerical tests in one dimension to show that the algorithm gives very good results for one 
dimensional flows in both adiabatic and isothermal MHD. For the one dimensional tests the 
field structure is either two or three dimensional. 

The algorithm has many astrophysical applications and is particularly suited to star for- 
mation problems. 
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1 INTRODUCTION 

Star forming regions are known to contain magnetic fields which 
are sufficiently strong to play an important part in the formation 
of the dense concentrations of matter which lead to stars. In order 
to describe the dynamics of such a system it is customary to begin 
with the equations of ideal magnetohydrodynamics (MHD). How- 
ever, the simplicity of the MHD equations is deceptive and hides the 
fact that there are numerous technical difficulties involved in their 
solution. Our aim in this paper is to describe a set of Smoothed Par- 
ticle Hydrodynamics (SPH, for a review see Monaghan 1992) equa- 
tions which overcome these difficulties and can be used to simulate 
MHD phenomena. We call the resulting algorithm Smoothed Parti- 
cle Magnetohydrodynamics (SPMHD). The equations appear to be 
accurate, robust and easy to apply and they do not suffer from the 
instabilities that exist in other formulations of the SPMHD equa- 
tions. 

An early application of SPH to MHD problems was to static 
magnetic polytropes iGingold & M onaghan 1977) who found good 
agreement with perturbation calculations. Dynamical problems 
were considered by Phillips 1 1983b D and appl i ed to st a r for- 
mation problems iPhillip J 1 1 982L 1 1 9833. 1 19851 Il986ibt bend 
1 1984 IPhillips & MonagharJ Il985l) . In the latter it was shown 
that when the conservation form of the equations was used 
an instability developed which took the form of SPH particles 



clumping. SPH blast waves in a magnetic medium were stud - 
ied by IStellingwerf & PeterkirJ ll99Cl Il994h. lHabe et alJ ll99lh . 
iMurrav. Wadslev & Bond ll99ot) and lMac Low et alj <1999h used 
a form of the SPH equations where the magnetic fields were up- 
dated on a grid and interpolated to the SPH particles. 

iMeglickil < 1994 1 19951) and lMeglicki et alj ll995l) used a for- 
mulation of SPMHD that uses a non-conservative (J x B) force, 
which is always stable and guarantees that the magnetic force is 
exactly perp endicular to the magne t ic fiel d. This formalism was 
also used by iBvleveld & Pongracid j 19961) and more recently by 
Cer queira and de G ouveia dal Pino (2001 and references therein) 
and H oskind 120021) . however the non-conservation of momentum 
leads to poor performance on shock-type problems. A conserva- 
tive fo rm of SPMHD h as been used bv | D olag 1 Bartelmann & Lesch| 
1 1999) and by Mari nho. And^eazza^^L^pinel ~ l200ll) since the 
magnetic field in their simulations remained in the regime where 
the instability does not appear. iMorriJ 1 199(3) suggested using a 
compromise between the conservative (tensor) force and the J x B 
formalism. Non-ideal MHD terms in SPH were also considered by 
iMorrisI ll996J) . who suggested using re sistive terms to control the 
divergence of the magnetic field and by Hosking (2002), who con- 
sidered the effects of ambipolar diffusion via a two-fluid model. 

The first technical difficulty with MHD simulations is that 
th e magnetic field comes w ith the constraint that V ■ B = 
0. IBrackbill & Barnes! Il980l) showed that in some finite dif- 
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ference codes, the failure to satisfy this constraint would lead 
to an instability. A number of different techniques have there- 
fore been developed to ensure that this constraint is satisfied. 
The first of these is to work with the vector potential A where 
B = V x A rather than wit h B. This approach was used 
by iGineold & Monaghar] 1 19771) for SPH si mulations. Others 
lEvans & Hawlevlll988UStone & Normanll99 j) construct their fi- 
nite difference equations so that, to the accuracy of the resolution, 
V ■ B = 0. Another approach commonly used is to clean up the 
magnetic field at every step by adding a gradient term to the com- 
puted field to produce a new field which satisfies the constraint. By 
giving up the conservation of momentum and including V ■ B terms 
in the momentum equation Powell 1 19941) (see lPowell et alll999l) 
produced a stable finite difference scheme for MHD (the 8 wave 
theory) which, however, appears to be less accurate for shocks. 
A comprehensi ve dis cussion of these and other schemes has been 
given by Toth 1 2000) who notes that even if V • B = most of 
these schemes produce magnetic forces which are not perpendic- 
ular to the exact m agn etic force J x B . Within the framework of 
SPH lBctrvd <200ll) (see lB0rve. Omang & Trulserl200ll) developed 
a non-conservative form of the SPMHD equations which have good 
stability properties and by the creation of closely spaced particles 
in regions with large spatial gradients it gives excellent a ccuracy. 

An alternative to any of these approaches is that of ljanhunerl 

(2000) who starts from the premise that non zero V ■ B terms 
may be generated but, if they are treated consistently, no instabili- 
ties wi l l occu r. The resulting set of equations has been derived by 
IPellaJ 1200 ll) starting with the rela tivistic formulatio n of gas dy- 
namics with electromagnetic fields. |janhunerJ (2000) showed that 
this formulation of MHD gas dynamics could be simulated using 
the HLL method ( Harten, Lax & Van Leer 1983) and he showed 
by thousands of test cases that positivity could be expected even if 
not proven. 

Our approach is to follow Janhunen 1 2000) and use equations 
which are consistent even if V ■ B does not vanish. To simulate 
shocks we introduce an artificial dissipation which guarantees that 
changes to entropy and thermal energy from viscous and ohmic 
dissipation are positive. The resulting set of equations conserves 
momentum and energy exactly. 

Another technical difficulty peculiar to SPH is that when 
a conservative force is used the SPH particles tend to clump 
in pairs in the presence o f tension. This was first noticed by 
Philli ps"& Monaehanl J 19851) and re-discovered by researchers ap- 
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2000). Several remedies have been proposed (e.g. 


Dvka. Randies & Ingel 


1997; 


Bonet & Kulasegaram 2000, 


2001) 



but they all either involve a significant increase in computation 
or cannot be applied where the particle configuration changes sig- 
nificantly. A remedy for the tensile instability which can be eas- 
ily app lied to astrophysical problems was proposed bv lMonaghanl 
(2000). The idea is to add a small artificial stress which prevents 
particles from clumping in the presence of a negative stress. This 
term has been shown to work well in elastic dynamics simulations 
iGrav et aljibol and we apply it here to the MHD case. We find 
that such a term very effectively removes the instability with few 
side effects. 

In f|2]we give the continuum form of the equations, and in Sj3] 
the SPH form of these equations. We construct appropriate dissi- 
pation terms for MHD in Sj4] The instability correction is discussed 
in Sj5]with details in the appendix. The time-stepping strategy is 
described in S|6| In [0we present the results of extensive numer- 
ical tests for one dimensional problems involving discontinuous 



initial conditions. In a companion paper JPrice & Monagharl2003l 
hereafter paper II) we derive the SPMHD equations from a varia- 
tional principle, including the case where the smoothing length is 
regarded as a function of local particle density. A self-consistent 
derivation of the SPMHD equations in the latter case is shown to 
increase the accuracy of SPMHD wave propagation. Two and three 
dimensional tests will be presented elsewhere. 



2 THE CONTINUUM EQUATIONS 

In the absence of dissipation the i th component of the acceleration 
equation is 



dv l _ 1 OS 13 
dt p dxi 



(1) 



where d/dt denotes the derivative following the motion, and the 
stress S 13 in the case of ideal MHD is defined by 



S v = _p 5 H + -L( B l B 3 - U l3 B 2 ), 
Mo 2 



(2) 



Here B 1 is the i th component of the magnetic field and po is the 
permittivity of free space. For SI units po — (4ty) / 10 7 . 

The time change of the magnetic field i s given by the induction 
equation. We follow Janhunen (2000) and Dellar (2001) and con- 
struct the induction equation so that it is consistent even if V • B 
does not vanish. The induction equation including ohmic dissipa- 
tion then becomes 



^ + V x (v x B) = -V x (r)J) - v(V ■ B), 



(3) 



where the last term is the magnetic current l Janhu nenl200bl FPellar 
2001 ) and J is the normal current 



J = fioV x B. (4) 

and rj is the magnetic diffusivity l/(apo) where a is the conduc- 
tivity. 

This induction equation can be written 

dB 
~~dt 



= (B ■ V)v - B(V ■ v) - V x (r}J). 



(5) 



This last form of the induction equation is the standard form when 
the constraint V • B = is used. Magnetic monopoles associated 
with V • B 7^ do not affect this equation. Taking the divergence 
of we find that monopoles evolve according to 



d_ 

dt 



(V • B) + V ■ (vV • B) 



0, 



(6) 



which has the same form as the continuity equation for the density 
and therefore implies that the volume integral of V - B is conserved. 

It is common to solve the acceleration equation with the ther- 
mal energy and continuity equations, but in this section we will 
assume that the thermal energy equation is replaced by the to- 
tal energy equation. Our aim is to derive a set of SPH equations 
which conserve total energy and momentum while ensuring that 
the change in entropy due to dissipation is positive. 

The total energy e per unit volume is defined by 



1 2 B 

e = P \ t; v + u H 

2 pno 



(7) 



where u is the thermal energy/unit mass. The equation for the rate 
of change of e can be written in terms of the stress according to 
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dt 



or as 

de 



(ev) + 



-V- 



B 2 



e + P + 



+ V-(Bx ( V J)) : 
B(v B) 



(8) ^ 



B x 



(9) 



2/xo y mo 

To derive SPH equations it is convenient to replace the en- 
ergy/unit volume by the energy per unit mass <T 

e = ~v 2 + u+ (10) 
2 2p^ 

The equation for e"is 

f=i^ + iv.(Bx (?? J)). (11) 

The thermal energy equation can be derived either from 1 101 
giving 



du 
IE 



de 
dt 



dv 



d 

dt 



B 2 



(12) 



dt dt \ 2/iop 

or by using the first law of thermodynamics including the ohmic 
heating term. Either way we find 

= v- v + V j\ 

dt p 

The final equation is the density equation 



(13) 



dt H 



(14) 



In addition we need an equation of state and an expression for 
the conductivity. In this paper we assume the gas is an ideal gas 
and we introduce an artificial dissipation into the SPH equations 
which we then interpret in terms of an artificial viscosity, thermal 
conductivity and magnetic diffusivity. 



3 THE SPH EQUATIONS 

We take as our fundamental equations the acceleration equation, 
the total energy equation and the density equation. To construct 
our SPMHD equations we follow the procedure used for the spe- 
cial relativistic equations IChow & Monaghar]|l997l) . Initially the 
equations will be set up assuming there is no ohmic dissipation. 
We will introduce artificial dissipation in the SPMHD equations in 
order to handle shocks and we will then show that to guarantee that 
the ohmic dissipation is always positive we must include an extra 
term in the induction equation. This extra term gives the appropri- 
ate extension of the induction equation to include the effects of an 
artificial conductivity. 

The acceleration equation for SPH particle a is (Monaghan 

1992) 



(15) 



— — = > m b I — =- H T + Ll ab ■-, 

dt ^ \pl p\ J dxi 

where W a b = W^(|r a — r b \, h) is the smoothing kernel. The dis- 
sipation term Tl ab will be discussed shortly. We write the energy 
equation in the absence of ohmic dissipation in the form 



dt~ V dxi \ p ) + p 2 dxi ' 



(16) 



which is similar t o tha t used for special relativistic SPH 
<Chow & Monagharj[l997il . The SPH equivalent of this equation 



\- (vis? vis* MO \ ow ab 

) m b \ = 1 = h il ab — — — , 



(17) 



where £l ab is a dissipation term analogous to H ab - Because of the 
symmetry of the terms in the summation total linear momentum 

m a v a , and energy ^2 m^e a are conserved. 

The induction equation in the absence of the ohmic term can 
be written 



' v ' \ dxi J p p \ Qx3 



dB l B J d 

The SPH form of this equation is 

dBl, 



, dp 

dxi 



(19) 



where vi denotes 

ba 



vi) . Equivalently we can use 



d Bl\ 1 ^ 4 R3 dW ab 

T7 — = — > m b v ba Bi——. (20) 
dt \p a I pi*-? dxi 

As is usual practice in SPH, the density is estimated via a sum- 
mation over neighbouring particles according to 



Pa 



b 



m b W ab , 



(21) 



or alternatively using the time derivative of this expression which 
gives the SPH form of the continuity equation 

dpq 

dt 



\ - j dW ab 
m b v ba - 



b 



dx' a 



(22) 



We note that equations I15i and 1 17t can be derived from a 
variational principle using 1 191 and 1221 as constraints, demonstrat- 
ing that these are indeed a consistent set of equations. This is pre- 
sented in paper II. 

The smoothing kernel we use is the usual spline-based kernel, 
given by 



W(q) 



a 

1? 



h 2 + h 3 , 



1(2 -q? 







0< q< 1; 
1 < q < 2; 

q > 2 



(23) 



where q = \r a — r b \/h, v is the number of spatial dimensions and 
the normalisation constant a is given by 2/3, 10/(77r) and 1/7T 
in 1, 2 and 3 dimensions respectively. The smoothing length h of 
particle a is set according to the usual rule 

h a oc f —J . (24) 

We implement this by evolving the smoothing length according to 
teenJl99dlMonagharll993) 

dh a _ _ha_d£a_ , 

dt up a dt ' 

This works extremely well for the tests presented in this paper since 
the density is evolved using the continuity equation I22i . We note, 
however, that the dependence of the smoothing length on the den- 
sity given by 1241 can be used to derive (again via a variational prin- 
ciple) a set of discrete equations for both SPH and SPMHD which 
self-consistently account for the extra terms which arise from this 
dependence. This set of equations is derived and implemented in 
paper II, where we demonstrate that it leads to increased accuracy 
in the propagation of MHD waves. In particular the formalism de- 
rived in paper II is natural to use when the density is calculated via 
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the SPH summation <2 It . In this paper we simply take the aver- 
age of the kernel to mainta in the symmetry in the mom entum and 
energy equations iHernaui st & KatJ ll989: Monaghan i 19921) . that 



Wab 



2 [W ab (h a ) 



+ W ab {h h )\ , 



and correspondingly 



dW ab 

dx l a 



dW ab {h a ) dWabQlb) 



dx a 



dx a 



(26) 



(27) 



4 DISSIPATION TERMS 

IChow & Monaghani 1 19971) discuss the SPH dissipation analo- 
gously to that associated with Riemann solvers. The key point is 
that the dissipation involves jumps in appropriate variables (mo- 
mentum, energy and density) between the left and right Riemann 
states multiplied by eigenvalues which can be interpreted as sig- 
nal velocities. In the SPH case we construct dissipation terms in a 
similar way. Thus 



n a6 



Kv si g(v a - v&) ■ j 

Pah 



(28) 



where K ~ 0.5 is a constant, v 3 i g is a signal velocity, and j is a unit 
vector from particle b to particle a. The density p ab = i(p a + pi) 
is an average density. 



J = 



Tab 

\r a b\ 



(29) 



In the relativistic case it was necessary to replace the velocity by 
the momentum and in the relativistic momentum use the velocity 
along the line of sight of the two particles in order to guarantee that 
the viscous dissipation makes a positive contribution to the thermal 
energy and therefore to the entropy. We will see that similar ideas 
are required here. 

The dissipation in the energy equation can be taken as 



Kv Bi g{e* a - e|)j 

Pab 



(30) 



where e* is an energy quantity which is related to e! Its precise form 
will now be deduced by considering the rate of change of thermal 
energy. 

From 1 1 21 using the SPH equations for the rate of change of 
velocity, energy, magnetic field and density we find that the mag- 
netic terms cancel leaving 

du a P a ^ i dWab 

dt ~ 



+ 



E 



2 A^ mbVab -f> T i 



Kv s 



m b - 



pab 



[(e*a - e*b) 

-(Va ' j)(Vo - Vb) ■ })]\ r ab\F ab , 



(3D 

where F a br a b ~ VWat, and F a b < is an even function of r a b- 

The first term on the right hand side is the adiabatic change 
in thermal energy due to the expansion or compression of the gas. 
This term does not change the entropy. The second term is the con- 
tribution to the change in thermal energy due to viscous dissipation, 
thermal conduction and ohmic heating. Only the first and last of 
these must contribute a non negative quantity to the change in the 
thermal energy. Heat conduction can either increase or decrease the 



thermal energy of an element of fluid. However, all three must con- 
tribute a positive quantity to the change of entropy of the system. 
The proof that the contribution to the entropy is positive is given in 
appendix[H] The terms due to viscous dissipation and ohmic dissi- 
pation must be negative definite because Fab < 0. 

It is natural to try and construct e* using the terms in?, namely 



- 1 2 B 

e = -v + u + . 

2 2/i p 



4.1 Viscous dissipation 



(32) 



The kinetic energy combined with the velocity terms in 1321 is not 
negative definite and in this form cannot be the correct viscous dis- 
sipation. We get a positive definite viscous dissipation in the veloc- 
ity terms by choosing 



e*a-el = ^(Va-i) 2 -^(Vb-i) 2 +Ua-U b +-!^—--^— 

2 2 2poPa 2/j.oPb 

We can then write the velocity terms in J3 li as 
1 



(33) 



l v 6 -JJ 



(34) 



which is negative definite. Recalling that F a b < the viscous con- 
tribution to the thermal energy is therefore positive. Note that the 
combination (e* — el) is not a simple difference because both e* 
and el involve j which depends on both particles. 

4.2 Ohmic dissipation 

The magnetic term is wrong because it can be positive or negative. 
In addition, it depends on the total field, whereas in a shock we 
would expect it to involve only the component perpendicular to the 
shock. We therefore replace the magnetic energy term by using the 
component of the field perpendicular to the line of sight of the two 
particles a and b, then we have 



* * 1 / .\2 1 / .\2 | 

e a -e b = -(v a -j) --(vi-j) +u a -u b 

+ [B 2 a - (B„ • j) 2 - B 2 b + (B b ■ i) 2 ] . (35) 

2pop ab L J 

The magnetic term is still not negative definite. To make it 
negative definite we need to add a term 

1 



■ {(B a ■ j) [(B a ■ j) - (B b • j)] - B a • [B a - B 6 ]} . (36) 

P0pa b 

With this new term added the magnetic contribution to the thermal 
energy becomes 



77^ [B 2 b — (B ab • j) 2 ] , 

2PabP0 L 



(37) 



which is negative definite and, when combined with F a b, gives a 
positive contribution to the thermal energy change. 

The interpretation of the extra magnetic terms is quite simple : 
when currents are present, and the conductivity is finite, the induc- 
tion equation requires an extra term as in ^5). The contribution to 
the rate of change of thermal energy from the new term is 



/in ^—^ 



Kv s 



nib- 



p0 * ' Pab 

b 



[Bab - ](Bab ■ j)] TabFab 



(38) 



The term which must be added to the induction equation for 
consistency can be deduced by noting that the expression for the 
rate of change of thermal energy J12I has a magnetic term 



Smoothed Particle Magnetohydrodynamics 5 



d ( B 



B dB B dp 



(39) 



dt V 2p/^o / p/io rfi P 2 / i o dt 

Comparing the first term with I38i we find that the SPH induction 
equation requires a term 

dB a 

dt 



Pa} m b 



b 



Kit 

— 2^(j x (B ab x j)r ab F ab , (40) 

Pab 



and we expect that the continuum version of this term should be 
some approximation to 

-Vx(i)VxB), (41) 
which when r\ is constant is 

77 [V 2 B - V(V • B)] . (42) 

By replacing the summation in I40i by an integral, and ex- 
panding in a Taylor series about r a , and assuming that v S i g is con- 
stant, we find that 1401 is proportional to 



Kv s igh 



V 2 B - -|v(V ■ B) 



(43) 



which is similar to the exact equation with ohmic diffusivity 77 oc 

KVsigh. 



4.3 The signal velocity 

We refer the read er to a general discussion of sig nal velocities in 
iMonashanl J 19971) and lChow & Monaghar] jl997f) . The key point 
is that it is the relative speed of signals from moving observers at 
the positions of particles a and b when the signals are sent along 
the line of sight. If there are no magnetic fields a good estimate of 
this signal velocity is 



Vsig = C a + C b — f3v ab ■ j, 



(44) 



where c a denotes the speed of sound of particle a and (5 ~ 1. 
The signal velocity is larger when the particles are approaching 
each other and in practice, the effects of shocks can be included 
by choosing = 2 (however when the artificial dissipation switch 
discussed in i|4,4l is used we find it is better to set j3 — 1 due to the 
stronger source term). If there are magnetic fields then a variety of 
other waves are possible. The fastest wave in a static medium along 
the x axis has speed 



B 2 2B*c 

c 2 H h -— = + 

ppo JPPa 



c 2 + 



B 2 2B*c 



(45) 



p/JO y/PPO 

A natural generalization of 1441 for the case of magnetic fields is to 
take 



v a +v b - [3v ab ■ j, 



(46) 



where 



1 

Va= 2 



o . B a 2B a ■ jc a 
ci H 1 ; + \ c& + 



B 2 

■L->a 



2B Q 



■ JCa 



\/PaPQ 



pa PQ \/PaP0 V pa PQ 

with a similar equation for v b . 
4.4 Artificial dissipation switch 

A switch to reduce the artificial viscosity away from shocks is given 
bv lMorris & Mona ghanH 1 9971) . Using this switch together with the 
suggestions of Balsara 1 1995) in multi-dimensional simulations can 
virtually eliminate the problematic effects of using an artificial dis- 
sipation in SPH. 



The key idea is to regard the dissipation parameter K (c.f. 
eauation l28l as a particle property. This can then be evolved along 
with the fluid equations according to 



dKa 

dt 



K a - K„ 



(48) 



such that in the absence of sources 5, K decays to a value K„ 
over a timescale r. The timescale r is calculated according to 

h 



(49) 



where h is the particle's smoothing length, v sig is the maximum 
signal propagation speed at the particle location and C is a dimen- 
sionless parameter with value 0.1 < C < 0.2. We conservatively 
use C = 0.1 which means that the value of K decays to K m i n over 
~ 5 smoothing lengths. 

The source term 5 is chosen such that the artificial dissipation 
grows as the particle approaches a shock front. We use 



5 = max(— V • v, 0), 



(50) 



such that the dissipation grows in regions of strong compression. 
Following Morris & Monaghan 1 1997) where the ratio of specific 
heats 7 differs from 5/3 (but not for the isothermal case), we multi- 
ply 5 by a factor 



In 



5/3- 



5/3- 1 



/ 



In 



7 + 1 

7-1 



(51) 



Note that o ur source term is a factor of two times larger than 
the term used by Mor ris & Monagfianl 1 19971) since our dissipation 
parameter K is required to be of order K ~ 0.5 at a shock front, 
whilst the usual SPH artificial viscosity parameter a is of order 
unity. We prefer this stronger source term since it provides suffi- 
cient damping in the Sod 1 1978) hydrodynamic shock tube prob- 
lem and in the MHD shock tube tests we describe in this paper (ie. 
K ma x ~ 0.5 for these problems). 

In order to conserve momentum the average value K — 
0.5 (id + Kb) is used in equation 1281 or J30I . A lower limit of 
Kmin = 0.05 is used to preserve order away from shocks (note 
that this is an order of magnitude reduction from the usual value of 
K — 0.5 everywhere). 

The numerical tests in ^demonstrate that use of this limiter 
gives a significant reduction in dissipation away from shocks whilst 
preserving the shock-capturing ability of the code. 



5 INSTABILITY CORRECTION 

The tensile insta bility is corrected via the method proposed by 

iMonaehanl fcOOOl) . The idea is add a small term which prevents 

particles clumping under negative stress. The momentum equation 

i J 51 becomes 
,(4TT 



dv a 
dt 



b 

+ R 



nib 



BiBj 



+ 



S, 
P 

BiBj 



dWab 



(52) 



where R is a function which increases as the particle separation 
decreases, given by 



R 



e ( W ab \ n 
2p \ W(Ap) J ' 



(53) 
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where W is the SPH kernel and W(Ap) is the kernel evaluated at 
the average particle spacing. Further details of the derivation of this 
term are given in appendix^ For all the simulations presented here 
the particles are setup with h — 1.5Ap and therefore in 1531 we 
compute the kernel in the denominator using Ap/h = 1/1.5. We 
use e = 0.8 and n = 5 throughout this paper and in one dimension 
we apply the correction only in the x-direction. 

Note that where the total energy equation 1 1 7t is used, the 
source term 



dt 



src £ 



m b v l a R 



BiBj 



+ 



BiBj 



dWab 
dXi,a 



(54) 



is added for consistency. 

We show in [Qthat this correction term very effectively re- 
moves the tensile instability with few side effects. 



6 TIMESTEPPING 

We integrate the SPMHD equations using a simple midpoint 
predictor-corrector method. Quantities are predicted according to 



A 1/2 = A" 
where 

A = [x,v x ,v y ,v z ,p,e,By,B z ,h,K} T , 



At (dA\ 
~2~ \~dt) 



-1/2 



(55) 



(56) 



with the energy e interchangeable for the thermal energy it. Note 
that we evolve the smoothing length alongside the particle equa- 
tions as discussed in j3|and that the dissipation parameter K is 
also evolved in accordance with the switch discussed in £14.4-1 The 
density is also included since in this paper the continuity equation 
1221 is integrated rather than using the density summation 

The rates of change dA / dt of these quantities are then com- 
puted via the SPH summations using the predicted values A 1 / 2 . 
The corrector step is given by 



A* = A + — f — ) , 
2 V dt J 

and 

A 1 = 2A* - A . 

The timestep is determined by the Courant condition 

h a b 



dt c = Ccourinin 



sig ,ab 



(57) 



(58) 



(59) 



where h a b = 0.5(h a + hb), the signal velocity v s i g is given by 
equation 1461 . except that we use 



Vsig.ab =Va+V b + f]\v ab • j| 



(60) 



with /3 = 1 when v a (,-j > (ie. where the dissipation terms are not 
applied). The minimum in 1591 is taken over all particle interactions 
and in this paper we use C cotlr = 0.8. 

Although this condition is sufficient for all of the simulations 
described here, in general it is necessary to pose the additional con- 
straint from the forces 



It a 

la J 



1/2 



dtf = min 
where a a is the acceleration on particle a. 



(61) 



7 NUMERICAL TESTS IN ONE DIMENSION 

The numerical scheme described in this paper has been tested 
on a variety of one dimensional problems. In order to demon- 
strate that SPMHD gives good results on problems involving dis- 
continuities in the physical variables we present res ults of Stan 
dard problems used to test grid-base MHD codes ( e.g.|Stoneet_al 
I l992t iDai & Woodward 1 1994 iRvu & JonesHl995t lBalsarj|l99 ~ 
bai & Woodwardl 19981) . The advantages of SPMHD are the sim- 
plicity with which these results can be obtained and the complete 
absence of any numerical grid. Further tests (MHD waves) are 
given in paper II since they incorporate the use of the variable 
smoothing length terms. 



7.1 Implementation 

The particles are allowed to move in one dimension only, whilst 
the velocity and magnetic field are allowed to vary in three dimen- 
sions. We use equal mass particles such that density changes cor- 
respond to changes in particle spacing. Unless otherwise indicated 
in this paper we integrate the continuity equation 122k the momen- 
tum equation 1151 . the total energy equation J17> and the induc- 
tion equation 1191 . This is the most efficient implementation of the 
SPMHD equations since it does not require an extra pass over the 
particles to calculate the density via the summation <2H . Similar 
results to those shown here are also obtained when the thermal en- 
ergy equation is integrated instead of the total energy. Additionally 
we note that whilst evolving the flux per unit mass 1201 instead of 
the flux density i 191 does not exactly maintain V • B = in one 
dimension, the associated errors are small and hence we also find in 
this case that the results are similar. Unless otherwise indicated the 
tests presented here are all performed with the artificial dissipation 
switch discussed in i|4. 41 turned on with minimum dissipation pa- 
rameter Kmin = 0.05. This results in very little dissipation away 
from shock fronts. 



7.1.1 Scaling 

The magnetic field variable is scaled in units such that the constant 
Ho is unity and numerical quantities are dimensionless. Note that 
the magnetic flux density B has dimensions 



[B 



[mass] 



[time] [charge] ' 
whilst no has dimensions 
'mass][length] 



[charge] 2 



(62) 



(63) 



Choosing mass, length and time scales of unity and specifying 
fio = 1 therefore defines the unit of charge. Re-scaling of the mag- 
netic field variable to physical units requires multiplication of the 
code value by a constant 



B 



physical 



Ho [mass] 
[length] [time] 2 



1/2 



B 



numerical • 



(64) 



For example, in cgs units, with mass, length and time scales of unity 
the magnetic flux density in Gauss is given by 



B C9S = (4tt) i/2 B 

numerical • 



(65) 
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7.1.2 Initial conditions 

Integration of the continuity equation 1221 requir es some sm ooth- 
ing of the initial conditions and we follow | Monaehan < 19971) such 
that when an initial quantity A is discontinuous it is smoothed ac- 
cording to the rule 



A = 



A L + A R e 



:/d 



1 



x I d 



(66) 



where Al and A R are the uniform left and right states with respect 
to the origin and d is taken as half of the largest initial particle sep- 
aration at the interface (ie. the particle separation on the low den- 
sity side). Note however that we do not smooth the initial velocity 
profiles except in the rarefaction test. Where the initial density is 
smoothed the particles are spaced according to the rule 



■ Xa-1) 



2p R A R 



(67) 



where A_r is the particle spacing to the far right of the origin with 
density p R . Note that initial smoothing lengths are set according to 
the rule h oc 1/p and are therefore also smoothed. Where the total 
energy e is integrated we smooth the basic variables u and B and 
construct the total energy from J I ()i . There is some inconsistency 
in the smoothing in this case as it is not possible to self consistently 
construct smooth profiles of p, u, B y , B z and e with the above 
smoothing. This can cause a small glitch in the initial conditions 
when the total energy equation is used. 

Such smoothing of the initial conditions can be avoided al- 
together if the density summation <2H is used, particularly if the 
smoothing length is updated self-consistently with the density. This 
is demonstrated in paper II. 



7.1.3 Boundaries 

Boundary conditions are implemented in one dimension by sim- 
ply fixing the properties of the 6 particles closest to each boundary. 
Where the initial velocities of these particles are non-zero their po- 
sitions are evolved accordingly and a particle is removed from the 
domain once it has crossed the boundary. Where the distance be- 
tween the closest particle and the boundary is more than the initial 
particle spacing a new particle is introduced to the domain. Hence 
for inflow or outflow boundary conditions the resolution changes 
throughout the simulation. 



7.1.4 Dissipative terms 

Direct application of the dissipation terms described in Sj4]provides 
no smoothing for the velocity in the y (and z) directions since we 
use (v a — vt) • j in equation 1281 and the particles are restricted to 
move along the x axis only. If the particles were allowed to move 
in the y(z) direction (with velocity v y (v z ) ) such smoothing would 
naturally be present. Therefore in the simulations presented in this 
paper we use 



Kv ai g(Va - V b ) ■ V 
pab 



(68) 



where v is a unit vector along the direction of the particle velocities 
given by 



v a - v b 

v = 1 T- 

|v„ - v&| 

The term in the momentum equation is then 



(69) 



b 



m b U ab \G ab 



(70) 



where we note that the gradient of the kernel may be written as 
VaWab = jGab, which we replace with vG ab where G ab is a 
scalar function (note that G ab — F ab h/\r ab \). 
In the dissipative energy we use 



e *a = Jj(«a ■ V) 2 +U a + i(B a ■ j) 2 , 



(71) 



with the contribution to the thermal energy equation from the ki- 
netic terms given by 

1 



[(v a • v) - (v 6 ■ v)]' 



(72) 



Note that the dissipative terms are only applied where v ab ■ j > 0. 
7.2 Simple advection test 

This simple test is described in lEvans & Hawlevl ll988h and in 
IStoneetaljil992l) and measures the ability of an algorithm to ad- 
vert contact discontinuities. A square pulse of transverse magnetic 
field is setup and adverted a distance of five times its width with 
the pressure terms switched off. The current density J is calculated 
in order to ascertain that the method does not produce sign rever- 
sals or anomalous extrema in this quantity. In SPH we compute this 
quantity using 



J a = V x B a = m b (B a - B b ) x V a W ab . 



(73) 



We perform this test simply by using a magnetic pressure that 
is negligible compared to the gas pressure. We setup 100 particles 
placed evenly along the x axis with constant velocity in the posi- 
tive x-direction and use a pulse that is initially 50 particle spacings 
wide. The pulse is not initially smoothed in any way and periodic 
boundary conditions are enforced using ghost particles (this is also 
a good test of the periodic boundary conditions since the particles 
are continually crossing the domain). 

The SPMHD results are shown in FigureQafter adverting the 
pulse a distance of five times its width (in this case 5 crossings of 
the computational domain). The top panel shows the results with 
the artificial dissipation terms turned off. The spread in the dis- 
continuities are kept to less than a particle spacing, showing no 
visible dispersion or diffusion whatsoever, suggesting that SPH in- 
deed handles contact discontinuities very well. The current density, 
which is analytically given by a delta function at each disconti- 
nuity, is also com puted very well by the SPH approximation (see 
lMonaghanlll992h . When the dissipative terms are turned on using 
the switch ( i)4.4t a small smoothing of the field is observed (bottom 
panels), however this still c ompares favourably with the schemes 
shown in St one et alJ l!992). 

7.3 Shock tubes 

The first s hock tube test we perform was first described by 
iBrio & Wii Jl988t) and is the MHD analog of the|SocJ < 19781) shock 
tube problem. The problem consists of a discontinuity in pressure, 
density, transverse magnetic field and internal energy initially lo- 
cated at the origin. As time develops complex shock structures 
develop which only occur in MHD because of the different wave 
types. Specifically the lBrio & WuUl988h problem contains a com- 
pound wave consisting of a slow shock attached to a rarefaction 
wave. The existence of such intermediate shocks was contrary to 
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[t] 



Figure 1. Results of the advection of a square pulse of transverse magnetic field 50 particle separations wide a distance of five times its width. In the absence 
of dissipative terms the discontinuities are kept to less than a particle spacing (top) and the current density (top right) shows almost no spread (analytically 
this is a delta function at each discontinuity). With the dissipative terms included with the Morris & Monaghan 1 1997) switch a small amount of smoothing is 
observed (bottom panels) 



the expectations of earlier theoretical studies iBrio & Wil l 19881. 
This problem is now a standard test for any astrophysical MHD 
code and has been used by many authors (e 1 gjStone j eyuj fl992t 
iDai & Woodwaro1l994fRvu & JoneJl995tlBalsariil998h 

We set up the problem using approximately 800 equal mass 
particles in the domain x = [—0.5,0.5]. Initial conditions to 
the left of the discontinuity (hereafter the left state) are given 
by (p, P,v x ,v y , By) = [1, 1, 0, 0, 1] and conditions to the right 
(hereafter the right state) are given by (p, P, v x , v y , B y ) = 
[0.125,0.1,0,0,-1] with B x = 0.75 and 7 = 2.0. The results 
are shown in Figure E| at time t = .1 and compare well with the 
numerical solution in Balsara 1 1998) (solid lines). Similar results 
to Figure|2|are obtained when the thermal energy equation is inte- 
grated. 

In the second shock tube test (Figure[3j, we demonstrate the 
usefulness of the artificial dissipation switch by considering a prob- 
lem which involves both a fast and slow shock. We consider the 
Riemann problem with left state (p, P, v x ,v y , B y ) = [1, 1, 0, 0, 1] 
and the right state (p, P, v x , v y , B y ) = [0.2 0.1, 0, 0, 0] with 
B x = 1 and 7 = 5/3 . This test h as been used bylDai & Woodward 
ll994) . lRvu & Jones! il99l and lBalsard ll998» and illustrates the 
formation of a switch-on fast shock. Similarly to the previous test 
we set up the simulation using approximately 800 particles in the 
domain x — [—0.5,0.5]. The results are shown in Figure [3] at 
time t = 0.15 and compare well with the exact solution given by 
lRvu& JonesNl995l) (solid lines). The advantages of the dissipation 



switch are apparent in this problem since it contains both a fast and 
slow shock. In a run with dissipation parameter K — 0.5 every- 
where the fast shock is significantly damped. In Figure [5] we see 
that the fast shock is resolved, although some small oscillations are 
observed behind the shock front. These oscillations can be removed 
entirely by using a slightly higher minimum dissipation parameter 
(Kmin = 0.1). 

The third test illustrates the formation of seven 
discontinuities in the same problem (Figure |4}- The 
left state is given by (p, P, v x , v y , v z , B y , B z ) = 
[1.08, 0.95, 1.2, 0.01, 0.5, 3.6/(4tt) 1/2 , 2/(4tt) 1/2 ] 
and the right state (p, P, v x , v y , v z , B y , B z ) = 
[l,l,0,0,0,4/(47r) 1/2 ,2/(47r) 1/2 ] with B x = 2/(4n) 1/2 
and 7 = 5/3. Since the velocity in the x-direction is non-zero at 
the boundary, we continually inject particles into the left half of the 
domain with the appropriate left state properties. The resolution 
therefore varies from an initial 700 particles to 875 particles at 
t = 0.2. The results are shown in Figure |4| at time t = 0.2. 
The SPMHD solution compares e xtreme ly well with the exact 
solution taken from iRvu & Joned ^995) (solid line) and may 
al so be com pared with the numerical solution in that paper and 
in iBalsaral ll998t) . The thermal energy and density profiles are 
slightly improved by our use of the total energy equation. Note 
that the initial velocity profiles are not smoothed for this problem, 
resulting in the small starting error at the contact discontinuity. 

The fourth test (Figure [5} is similar to the previ- 
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x x 

Figure 2. Results of the lBrio & Wulll988t) shock tube test. To the left of the origin the initial state is (p, P,v x ,v y ,B y ) = [1, 1, 0, 0, 1] whilst to the right 
the initial state is (p, P, v x ,v y , B y ) = [0.125, 0.1, 0, 0, —1] with B x = 0.75 everywhere and 7 = 2.0. Prof iles of density, p ressure, v x , v y , thermal energy 
and By are shown at time t = 0.1. Points indicate the SPMHD particles whilst the numerical solution from Balsara 1 1998) is given by the solid lines. The 
artificial dissipation switch with K m i n = 0.05 is used. 
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Figure 3. Results oftheMHD shock tube test with left state (p, P, v x , v y , B y ) = [1,1,0,0,1] and the right state (p, P,v x ,v y , B y ) = [0.2,0.1,0,0,0] 
with B x = 1 and 7 = 5/3 at time t = 0.15. The problem illustrates the formation of a switch-on fast shock and the solution contains both a fast and slow 
shock. Solid points indicate the SPMHD particles whilst the exact solution is given by the solid line. The artificial dissipation switch is used. Without this 
switch the fast shock is significantly damped. 
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Figure 4. Results of the MHD shock tube test with left state (p, P, v x , v v , v z , B y , B z ) = [1.08,0.95, 1.2, 0.01,0.5, 3.6/(4tt) 1 / 2 , 2/(4vr) 1 / 2 ] and right 
state (p,P,v x ,v y , v z , B y , B z ) = [1, 1, 0, 0, 0, 4/(4tt) 1 / 2 , 2/(4tt) 1 / 2 ] with B x = 2/(4tt) 1 /2 an d 7 = 5/3 at time i = 0.2. This problem illustrates 
the formation of seven discontinuities. The exact solution is given by the solid line whilst points indicate the positions of the SPMHD particles. The artificial 
dissipation switch is used with K min =0.05 
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Figure 5. Results of the isothermal MHD shock tube test with left state (p, v x 



state (p, P, i) x 



,B y ,B z ) = [1.08, 1.2, 0.01, 0.5, 3.6/(4tt) 1 / 2 , 2/(4tt) 1 / 2 ] and right 



, B y ,B z ) = [1, 0, 0, 0, 4/(4vr) 1 /2 i 2/(4tt) 1 /2] wml B x = 2/(4w) 1 / 2 and an isothermal sound speed of unity at time t = 0.2. This 



problem illustrates the formation of six discontinuities in isothermal MHD. Solid points indicate the position of the SPMHD particles which may be compared 
with the exact solution given by the solid line. The artificial dissipation switch is used with K m i n = 0.05. 
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Figure 6. Results of the MHD shock tube test with left state (p, P, v x , v y , B y ) = [1, 1, -1, 0, 1] and right state (p, P, v x , v y , B y ) = [1,1,1,0,1] with 
B x = and 7 = 5/3 at time t = 0.1. This problem illustrates the formation of two magnetosonic rarefactions. The exact solution is given by the solid line 
whilst points indicate the position of the SPMHD particles. The artificial dissipation switch is used with K m i n = 0.05. 



ous version except that an isothermal equation of state is 
used. The left state is given by (p, v x , v y , v z , B y , B z ) = 
[1.08, 1.2, 0.01, 0.5,3. 6/(47r) 1/2 ,2/(47r) 1/2 ] and the right state 
(p,v x ,v y ,v z ,B y ,B z ) = [l,0,0,0,4/(47r) 1/2 ,2/(47r) 1/2 ] with 
B x = 2/(47r) 1//2 and an isothermal sound speed of unity. Results 



are shown in Figure|5|at time t = 0.2 and co mpare very well with 
the numerical results given in lBalsardil993) (solid line). 

The fifth test shows the formation of two magnetosonic 
rarefactions. The left state is given by (p, P, v x , v y , B y ) = 
[1,1,-1,0,1] and the right state by (p, P, v x , v y , B y ) = 



14 Price & Monaghan 



[t] 



t=0.03 



2000 



1500 



1000 - 



500 



-0.5 



0.5 



-0.5 



0.5 



40 



20 



-20 - 



-40 



- 



-0.1 - 



-0.2 



-0.5 



0.5 




0.02 - 



^0.02 



-0.04 = 



-0.06 



-0.5 



0.5 




1 .2 



600 - 



0.8 



400 



0.6 



0.4 



200 



-0.5 



0.5 



-0.5 



0.5 



Figure 7. Results of the MHD shock tube test with left state (p, P,v x ,v y ,v z ,B y ,B z ) = [1, 1, 36.87, -0.155, -0.0386, 4/(4tt) 1 '' 2 , 1/(4tt) 1 / 2 ] and right 
state (p, P,v x ,v y ,v z ,By,B z ) = [1, 1, -36.87, 0, 0, 4/(4tt) 1 / 2 , 1/(4tt) 1 / 2 ] with B x = 4.0/(4tt) 1 / 2 and 7 = 5/3. Results are shown at time t = 0.03. 
This problem illustrates the formation of two extremely strong fast shocks of Mach number 25.5 each. Solid points indicate the position of the SPH particles 
whilst the exact solution is given by the solid line. The artificial dissipation switch is used with K m i n = 0.05. The overshoots in density, pressure and 
magnetic field are a result of our integration of the continuity equation and neglect of terms relating to the gradient of the smoothing length, (these terms are 
derived in paper II). 
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[1, 1, 1, 0, 1] with B x = and 7 = 5/3. Results are shown in 
Figure |6| at time t = 0.1 and compare extremely well with the ex- 
act solution from Rvu & Jones! (1995) (solid line). Outflow bound- 
ary conditions are used such that the resolution varies from an 
initial 500 particles down to 402 particles at t = 0.1 in the do- 
main x = [—0.5,0.5]. The artificial dissipation switch is turned 
on with Kmin = 0.05 although very little dissipation occurs in 
this simulation since the artificial dissipation is only applied for 
particles approaching each other. With unsmoothed initial condi- 
tions we therefore observe some oscillations behind the rarefaction 
waves, which are removed in t his case by s mooth ing the initial dis- 
continuity slightly. As noted in Monaghan 1 1997) use of the density 
summation also improves the results for this type of problem. 

The fin al tes t, tak en from iDai & Woodwar d 

Jl994l) and iBalsad 1 1998), illustrates the formation 
of two fast shocks, each with Mach number 25.5. 
The left state is given by (p, P,v x ,v y ,v z , B y , B z ) = 
[1, 1, 36.87, -0.155, -0.0386, 4/(4tt) 1/2 , 1/(4tt) 1/2 ] 
and the right state by (p, P, v x , v y , v z , B y , B z ) = 
[1,1, -36.87, 0,0,4/(47r) 1/2 ,l/(47r) 1/2 ] with B x 
4.0/(4tt) 1/2 and 7 = 5/3. Results are shown in Figure at 
time t = 0.03. Inflow boundary conditions are used such that 
the resolution varies from an initial 400 particles up to 1286 
particles at t — 0.03 in the domain x = [—0.5, 0.5]. The artificial 
dissipation switch is turned on with K m i n = 0.05. The results 
co mpare extremely w e ll with the exact solution (solid line) given 
by |DaT & Woodward ( 1994) and with the numerical solution 
given by IDai & Woodward 1 19941) and lBalsaral ll998l) . especially 
given the extreme nature of the problem. The spikes in transverse 
velocity components are starting errors due to the fact that for this 
problem we do not smooth the initially discontinuous velocity 
profiles in any way. There is some advantage to integrating the 
total energy equation for this type of problem since using the 
thermal energy equation produces a large spike in thermal energy 
at the discontinuity and numerical noise behind the shocks. The 
overshoots in density, pressure and magnetic field are a result of 
our integration of the continuity equation and neglect of terms 
relating to the gradient of the smoothing length. These terms are 
derived and implemented in paper II and are shown to remove the 
errors seen here. 



8 SUMMARY 

We have shown how SPH equations for MHD can be formulated 
with the following features: 



(i) The equations us e the continuum equations of Ijanhunerl 
( 2000) and Dellar 1 2001) which are consistent even when the diver- 
gence of the magnetic field is non zero. Consequently, even though 
non zero V ■ B may be produced during the simulation, it is treated 
consistently. We find that insisting on consistency with fundamen- 
tal principles is the key to deriving stable SPH equations. 

(ii) The equations contain artificial dissipation. We require these 
dissipation terms to result in positive definite changes to the en- 
tropy and this places strong constraints on the form of the dissipa- 
tion. Our equations guarantee that the resulting viscous and ohmic 
dissipation produces positive changes in the thermal energy . In 
ensuring these equations are consistent with the fundamental re- 
quirement that the entropy should increase we are led to introduce 
a term in the induction equation which is anal ogous to the induction 
equation for non ideal MHD. The use of the lMorris & Monaghanl 



( 1997) switch very effectively minimizes the effect of the artificial 
dissipation away from shocks. 

(iii) The SPMHD equations also incorporate a simple technique 
to prevent an instability due to the tension arising from the magnetic 
stress. 

The resulting equations, when implemented with a simple pre- 
dictor corrector scheme, give good results for a wide range of shock 
tube problems. While we have yet to apply our algorithm to prob- 
lems in two and three dimensions the present results encourage us 
to believe that our SPMHD code will provide a secure basis for 
astrophysical MHD problems. 



APPENDIX A: ARTIFICIAL STRESS 

In this appendix we describe some of the details of the artificial 
stress required to prevent clumping. We diagonalise the magnetic 
part of the stress tensor by rotating the coordinate system so that the 
z axis lies along the magnetic field. The magnetic field is then B' = 
(0, 0, B) and the stress tensor is has non zero components M' xx = 
-B 2 /(2 A t ), M' yy = -B 2 /(2 A t ) , and M' zz = B 2 /{2p ). The 
sign of the first two is associated with compression and the sign of 
the third is associated with tension. To remove the tension term at 
close range we add a term to M zz so that it is negative when the 
particles approach . The term we choose is RB 2 , where 



2p \W(A P ) J 



(Al) 



where e ~ 0.4 and n ~ 4 and in the tests shown in this paper we 
have e = 0.8 and n — 5. The precise values of these quantities is 
not important. W a b is the kernel and W(Ap) is W evaluated at the 
average particle spacing. 

Rotating back to the original coordinates we find that the term 
we subtracted is equivalent to defining a new magnetic stress 



ML = Ma + RBiBj 



(A2) 



APPENDIX B: POSITIVITY OF THE ENTROPY CHANGE 

In this appendix we demonstrate that the dissipation terms intro- 
duced in ij4|lead to a positive definite increase in the entropy. 

The second law of thermodynamics shows that the change of 
entropy per unit mass s a of particle a is given by 



ds a du a Pa dp a 



dt 



dt 



pi dt 



(Bl) 



where T a is the temperature (absolute) of particle a. 

From <31> . 1341 . and 1371 . and noting that the second term of 
JB1L when expressed in SPH form, cancels the first term of <3U . 
we find that 



dSq 

' dt 



b 



m b Kv sig 

Pab 



(-K ■ j - Vf, ■ j] 2 



(B2) 



— [B 2 ;, — (B ab ■ j) 2 ] + U a — U b J r ab F a 



POpab 

Because F a b < Owe can rewrite this equation as 



rp ds a _ \ - m b Kv slg 



Pa b 



(u a - u b )r ab F ab , 



(B3) 
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where Q a > is the contribution to the entropy change from the 
viscous and ohmic dissipation which we have shown is positive def- 
inite. The second term is a heat conduction term which can increase 
or decrease the thermal energy of a particle but it must not result 
in a decrease of the entropy of the system. Note that if u a > Ub 
then the thermal conduction causes a decrease in the thermal en- 
ergy of particle a. This form of the conduction term that a rises here 
is similar to that derived by IClearv & Monaghan! J 19991) for SPH 
heat conduction problems. 

The change in the total entropy S is then 



dS v-^ ds a 

Tt = E™°^r 

a 



Ca 



T a 

a 

where 

_ m b Kv alg r ab F ab 

Pab 



(u a — u b ) 

J- a 



< 0. 



(B4) 



(B5) 



We can interchange the labels in the second term of <B4> and 
combine with the original form to write this term as 



sEE^^-^^-fr 



(B6) 



If, as is normally the case, u is a monotonically increasing function 
of T then, if T a > T b we have u a > Ub and 



(B7) 



so that the second term on the right hand side of <B4> is positive. 
The change in the entropy due to thermal conduction is therefore 
positive. 
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